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ABSTRACT 


The  specific  objective  of  this  research  has  been  the  development  and  improvement 
of  theoretical  models  to  simulate  the  effect  of  nonsphericity  on  single-scattering  properties 
of  cirrus  cloud  particles  in  the  visible  and  infrared  spectral  regions.  First,  we  have  shown 
that  using  a  matrix  inversion  scheme  based  on  a  special  LU  factorization  rather  than  on  the 
standard  Gaussian  elimination  significantly  improves  the  numerical  stability  of  T-matrix 
computations  for  nonabsorbing  and  weakly  absorbing  nonspherical  particles.  Second,  we 
use  exact  T-matrix  computations  and  the  Kirchhoff  approximation  to  show  that  the  ^-function 
transmission  peak  predicted  by  the  GO  approximation  for  hexagonal  ice  crystals  is  an 
artifact  of  GO  completely  ignoring  physical  optics  effects  and  must  be  convolved  with  the 
Fraunhofer  pattern,  thereby  producing  a  phase  function  component  with  an  angular  profile 
similar  to  the  standard  diffraction  component.  Third,  we  have  used  the  improved  T-matrix 
method  to  compute  the  linear  depolarization  ratio  for  polydispersions  of  randomly  oriented 
ice  spheroids,  circular  cylinders,  and  Chebyshev  particles  with  sizes  typical  of  young 
contrails.  We  have  shown  that  ice  crystals  with  effective  radii  as  small  as  several  tenths  of 
a  micron  can  already  produce  6  exceeding  0.5  at  visible  wavelengths. 
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INTRODUCTION 


This  research  project  was  intended  to  deliver  crucial  modeling  products  for  use  in  the 
project  headed  by  Dr.  B.-C.  Gao  and  called  "Cirrus  Cloud  Characterization  and  Correction." 
It  is  well  known  that  the  accuracy  of  determining  the  cirrus  cloud  optical  thickness  and 
correcting  the  effect  of  cirrus  on  the  retrieved  aerosol  and  lower-cloud  optical  thicknesses 
and  surface  reflectance  critically  depends  on  the  accuracy  of  modeling  the  cirrus  cloud 
particle  single-scattering  properties,  first  of  all  the  scattering  phase  function.  The  knowledge 
of  the  phase  function  is  also  crucially  important  in  an  accurate  evaluation  of  the  radiative 
effect  of  cirrus  on  climate.  Cirrus  clouds  consist  of  ice  particles  having  an  wide  variety  of 
shapes.  For  typical  cirrus  environments  ice  crystals  exist  both  in  single  and  aggregated 
forms,  and  the  surface  of  the  ice  particles  can  be  highly  irregular.  Because  of  ignorance 
about  the  distribution  of  cirrus  particles  over  shapes,  knowledge  of  cirrus  scattering  phase 
function  has  progressed  very  slowly.  Theoretical  phase  function  calculations  have  mostly 
assumed  that  cirrus  ice  particles  are  spheres  or  randomly-oriented,  infinitely-long  circular 
cylinders;  this  remains  a  common  assumption  both  in  global  climate  models  and  in  satellite 
retrieval  algorithms.  The  geometric-optics  hexagonal-crystal  phase  functions  of  Takano  and 
Liou  have  become  quite  popular,  but  do  not  seem  to  fit  the  known  facts  very  well.  Phase 
function  calculations  for  small  ice  particles  have  been  difficult  and  computationally 
intensive. 

Because  of  these  factors,  the  focus  of  this  research  during  the  period  23  September 
1996  -  30  September  1 998  was  on  further  development  and  use  state-of-the-art  numerical 
techniques  for  computing  scattering  properties  of  nonspherical  ice  crystals.  The  specific 
objective  has  been  the  development  and  improvement  of  theoretical  models  to  simulate  the 
effect  of  nonsphericity  on  single-scattering  properties  of  cirrus  cloud  particles  in  the  visible 
and  infrared  spectral  regions. 
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IMPROVEMENTS  OF  THE  T-MATRIX  METHOD 

The  Waterman's  T-matrix  approach  is  one  of  the  most  powerful  exact  techniques  for 
computing  light  scattering  by  nonspherical  particles  based  on  solving  Maxwell's  equations. 
However,  even  this  method  can  exhibit  convergence  problems  when  any  of  the  variables 
defining  the  scattering  particle  (size  parameter,  deviation  from  sphericity,  refractive  index) 
becomes  too  extreme.  Standard  T-matrix  computations  become  especially  ill-conditioned 
for  particles  with  a  small  or  zero  imaginary  part  of  the  refractive  index  because  of  the  strong 
effect  of  the  ripple  structure.  For  example,  the  maximum  convergent  equivalent-sphere  size 
parameter  (i.e.,  the  ratio  of  particle  circumference  to  wavelength  of  scattered  light  for  the 
surface-equivalent  sphere)  in  double-precision  FORTRAN  computations  for  oblate  spheroids 
with  refractive  index  1.53  and  aspect  ratio  2  is  only  8,  whereas  the  maximum  convergent 
size  parameter  for  the  same  particles  but  with  a  moderately  absorbing  refractive  index  of 
1.53  +  0.001  i  is  33. 

The  origin  of  the  numerical  instability  of  the  standard  T-matrix  procedure  for  extreme 
values  of  particle  characteristics  can  be  explained  as  follows.  Calculations  based  on  the 
extended  boundary  condition  method  (EBCM)  assume  the  representation  of  the  T-matrix  in 

the  form  T - Q~^  RgQ,  where  the  elements  of  the  matrices  Q  and  RgQ  are  integrals  over 

the  particle  surface.  The  numerical  inversion  of  the  matrix  Q  is  usually  performed  using  the 
standard  Gaussian  elimination  (GE).  Unfortunately,  the  calculation  of  the  inverse  matrix 
Q~'  is  an  ill-conditioned  procedure  strongly  affected  by  round-off  errors  and  by  the  fact  that 
different  elements  of  the  Q  matrix  can  differ  by  many  orders  of  magnitude.  As  a  result,  T- 
matrix  computations  for  extreme  particle  parameters  can  be  poorly  convergent  and  even 
divergent. 

This  sensitivity  of  the  standard  T-matrix  procedure  to  weak  or  zero  absorption  has 
been  a  serious  limiting  factor  since  many  commonly  encountered  substances  are  weakly 
absorbing  in  some  spectral  ranges.  The  example  of  water  ice  in  the  visible  is  particularly 
important  because  quite  often  a  significant  fraction  of  cirrus  and  contrail  ice  particles  are  not 
much  larger  than  a  visible  wavelength,  thus  potentially  precluding  the  use  of  the  geometric 
optics  approximation  (GO)  in  light  scattering  computations. 
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During  the  course  of  this  research,  we  have  significantly  improved  the  T-matrix 
technique  for  computing  light  scattering  by  nonspherical  ice  crystals.  Specifically,  we  have 
shown  that  using  a  matrix  inversion  scheme  based  on  a  special  LU  factorization  rather  than 
on  the  standard  Gaussian  elimination  significantly  improves  the  numerical  stability  of  T- 
matrix  computations  for  nonabsorbing  and  weakly  absorbing  nonspherical  particles.  We 
have  also  developed  an  improved  scheme  for  evaluating  Clebsch-Gordon  coefficients  with 
large  quantum  numbers  which  allowed  us  to  extend  the  analytical  orientational  averaging 
method  developed  by  Mishchenko  [J-  Opt.  Soc.  Am.  A  8,  871  (1991)]  to  larger  size 
parameters.  As  a  result,  the  maximum  convergent  size  parameter  for  particles  with  small 
or  zero  absorption  can  increase  by  a  factor  of  several  and  can  exceed  100.  The  ability  of 
the  new  T-matrix  scheme  to  treat  such  large  size  parameters  is  accompanied  by  an  extremely 
high  numerical  efficiency  which  makes  our  code  orders  of  magnitude  faster  than  any 
alternative  technique  for  exactly  computing  light  scattering  by  nonspherical  particles  in 
random  orientation.  The  unique  capabilities  of  the  T-matrix  code  make  it  very  useful  in 
practice,  but  also  make  difficult  checks  of  its  numerical  accuracy  since  independent  results 
for  the  largest  size  parameters  cannot  be  obtained  with  any  other  currently  available 
method.  Therefore,  we  have  paid  special  attention  to  making  sure  that  our  calculations  fully 
satisfy  such  fundamental  physical  constraints  as  symmetry,  reciprocity,  and  energy 
conservation.  Also,  we  have  computed  benchmark  results  for  a  challenging  test  case  that 
can  be  used  for  checking  the  accuracy  of  the  most  advanced  nonspherical  scattering  codes 
at  higher  frequencies. 

Using  the  present  version  of  the  T-matrix  code,  we  have  extended  comparisons  of 
exact  T-matrix  and  approximate  ray  tracing  calculations  to  much  larger  size  parameters  and 
to  all  elements  of  the  scattering  matrix.  Our  results  suggest  that  equivalent-sphere  size 
parameters  larger  than  about  80  are  already  big  enough  to  ensure  acceptable  accuracy  of 
GO  phase  function  computations  (except,  perhaps,  at  exactly  the  backscattering  direction). 
However,  GO  calculations  of  the  other  elements  of  the  scattering  matrix  are  stronger 
affected  by  wave  effects  and  become  reasonably  accurate  only  at  significantly  larger  size 
parameters.  GO  calculations  of  lidar  depolarization  can  be  expected  to  be  especially 
inaccurate  unless  the  equivalent-sphere  size  parameter  exceeds  several  hundred. 
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An  interesting  result  of  our  calculations  is  that  the  T-matrix  method  can  be 
successfully  applied  to  large  sharp-edged  particles  such  as  finite  circular  cylinders  with 
equivalent-sphere  size  parameters  exceeding  100.  It  has  been  often  claimed  that  the 
presence  of  sharp  edges  can  be  difficult  to  handle  with  a  method  that  uses  "smooth" 
spherical  functions  in  the  internal  and  scattered  field  decompositions.  We  have  found, 
however,  that  the  use  of  a  special  numerical  integration  scheme  for  computing  the  surface 
integrals  needed  to  calculate  the  T  matrix  ameliorates  the  problem  of  sharp  edges  and  makes 
T-matrix  computations  for  cylinders  almost  as  accurate  as  those  for  surface-  and  aspect-ratio- 
equivalent  smooth-shaped  spheroids.  Much  more  difficult  problems  are  encountered  when 
the  T-matrix  method  is  applied  to  particles  with  large  aspect  ratios.  In  this  case  a  single 
spherical  function  expansion  of  the  internal  and  scattered  fields  can  fail,  and  the  use  of 
several  overlapping  subdomain  spherical  function  expansions  may  become  necessary. 

We  published  a  user  guide  describing  in  detail  a  software  implementation  of  the 
current  version  of  the  T-matrix  method  for  computing  light  scattering  by  polydisperse, 
randomly  oriented,  rotationally  symmetric  particles.  The  FORTRAN  T-matrix  codes  are 
publicly  available  on  the  World  Wide  Web  at  http://www.giss.nasa.gov/~crmim.  We 
provided  all  necessary  formulas,  described  input  and  output  parameters,  discussed  numerical 
aspects  of  T-matrix  computations,  demonstrated  the  capabilities  and  limitations  of  the  codes, 
and  discussed  the  performance  of  the  codes  in  comparison  with  other  available  numerical 
approaches. 

This  research  was  published  in  Refs.  1  and  2  below. 

INCORPORATION  OF  PHYSICAL  OPTICS  EFFECTS  INTO 
RAY-TRACING  PHASE  FUNCTIONS  COMPUTED  FOR  HEXAGONAL  ICE 

CRYSTALS 

It  is  well  known  that  a  convenient  way  of  representing  the  scattering  phase  function 
P(0)  for  aerosol  and  cloud  particles  is  expanding  it  in  Legendre  polynomials  as 
where  0  is  the  scattering  angle,  P„{cos  0)  are  Legendre  polynomials,  and  the  value  of  the 
upper  summation  limit  depends  on  the  desired  numerical  accuracy  of  the  expansion. 
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(1) 


max 

P(0)  =  X„  P„(COS0), 

n  =0 

Since  the  number  of  numerically  significant  terms  in  the  Legendre  expansion  is  finite  and 
often  relatively  small,  this  expansion  can  be  used  for  efficiently  computing  the  phase 
function  for  essentially  any  number  of  scattering  angles  with  a  small  consumption  of  CPU 
time.  Furthermore,  the  Legendre  expansion  coefficients  x„  can  be  used  to  directly  compute 
the  Fourier  components  of  the  phase  function  via  simple  and  exact  analytical  formulas, 
which  is  the  first  step  in  radiative  transfer  computations  using  different  numerical  techniques. 

The  Legendre  expansion  coefficients  for  the  widely  used  Henyey-Greenstein 
phase  function  are  given  by  the  simple  analytical  formula 

=  (2n  +  1)g",  (2) 

where  g  is  the  asymmetry  parameter.  Efficient  analytical  methods  based  on  solving 
Maxwell's  equations  exist  for  computing  the  expansion  coefficients  for  spherical  particles 
and  randomly  oriented,  rotationally  symmetric  nonspherical  particles.  For  irregular  particles 
with  sizes  much  larger  than  the  wavelength  of  the  incident  radiation,  such  as  cirrus  cloud 
particles  in  the  visible,  direct  numerical  solutions  of  Maxwell's  equations  do  not  currently 
exist.  Therefore  the  expansion  coefficients  have  to  be  computed  using  an  approximate 
technique  such  as  the  geometric  optics  approximation  (GO).  Using  the  orthogonality 
property  of  Legendre  polynomials,  we  easily  derive  from  equation  (1) 

rr 

d0P(0)Pn(cos0)sin0.  (3) 

The  integral  in  equation  (3)  can  be  calculated  numerically  by  using  a  quadrature  formula 
provided  that  the  phase  function  values  at  the  division  points  are  known.  This  numerical 
approach  works  well  if  the  phase  function  is  rather  smooth  but  becomes  problematic  for 


= 


2n  +  1 


i 
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particles  having  parallel  planes  such  as  hexagonal  columns  and  plates  or  finite  circular 
cylinders.  In  this  case,  the  standard  GO  predicts  a  strong,  infinitesimally  narrow  peak  in  the 
exact  forward-scattering  direction  which  is  caused  by  rays  that  undergo  two  refractions 
through  parallel  plane  facets  and  is  superimposed  on  the  diffraction  component  of  the  phase 
function.  This  effect  was  called  by  Takano  and  Liou  the  J-function  transmission. 

It  is  obvious,  however,  that  GO  predicts  the  infinitesimally  narrow  J-function 
transmission  peak  only  because  it  completely  ignores  physical  optics  effects.  Simple 
physical  optics  considerations  cause  us  to  conclude  that  although  a  strong  non-diffraction 
forward-scattering  peak  does  exist  and  can  be  qualitatively  explained  in  GO  terms  as  a 
manifestation  of  the  (5-function  transmission,  it  nonetheless  has  an  appreciable  angular  width 
comparable  to  that  of  the  Fraunhofer  diffraction  peak  and  a  diffraction-like  angular  profile. 
It  is  clear  that  however  large  a  particle  is  compared  to  the  wavelength,  physical  optics  effects 
will  preclude  the  appearance  of  perfect  singularities  in  the  scattering  pattern  like  the  6- 
function  transmission  peak  in  the  phase  function.  Instead,  a  wave  front  emerging  from  a  flat 
crystal  facet  should  spread  and  produce  an  angular  intensity  distribution  in  the  far-field  zone 
similar  to  the  well  known  Fraunhofer  diffraction  pattern.  Of  course,  in  the  theoretical  limit 
of  an  infinite  size  parameter  the  (^-function  transmission  peak  becomes  a  true  ^-function. 
However,  the  angular  width  of  the  ^-function  transmission  peak  is  always  comparable  to  that 
of  the  Kirchhoff  diffraction  component,  and  as  long  as  the  latter  is  computed  explicitly,  so 
should  be  the  tf-function  transmission  component. 

We  used  exact  T-matrix  computations  for  rather  large  nonspherical  particles  to 
demonstrate  that  the  effect  that  can  be  interpreted  in  geometric  optics  terms  as  (J-function 
transmission  through  parallel  planes  indeed  results  in  a  quasi-Fraunhofer  forward-scattering 
peak  rather  than  in  a  true  J-function  peak.  We  then  developed  a  very  simple  numerical 
procedure  which  incorporates  this  physical  optics  effect  in  the  standard  ray  tracing 
procedure  for  computing  the  phase  function  for  large  particles  having  parallel  plane  facets. 
This  procedure  not  only  makes  ray  tracing  computations  more  physically  relevant,  but  also 
simplifies  and  makes  more  accurate  the  computation  of  the  phase  function  and  its  Legendre 
expansion.  Although  the  accuracy  of  our  procedure  cannot  be  assessed  directly  due  to  the 
lack  of  exact  theoretical  methods  based  on  solving  Maxwell's  equations  and  applicable  to 
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size  parameters  exceeding  several  hundred,  our  approach  is  physically-based  and  appears 
to  be  simple  and  well  justified  since  it  consists  of  directly  computing  the  amount  of  energy 
contained  in  the  (J-function  transmission  peak  and  convolving  it  with  the  Fraunhofer  angular 
pattern. 

This  research  was  summarized  in  Ref.  3  below. 

Depolarization  of  lidar  returns  by  small  ice  crystals 

Because  of  intensifying  air  traffic,  there  has  been  increasing  interest  in  the  potential 
impact  of  aircraft  condensation  trails  (contrails)  on  climate  through  a  direct  radiative  forcing. 
The  estimation  of  the  climatic  effect  of  contrails  requires  knowledge  of  their  radiative 
properties,  which  can  be  substantially  different  from  those  of  ambient  cirrus  clouds.  The 
latter,  in  turn,  necessitates  the  determination  of  the  size  and  shape  distribution  of  contrail 
particles  and  their  time  evolution  using  in  situ  and/or  remote  sensing  techniques. 

Measurements  of  the  lidar  linear  depolarization  ratio  6  can  be  a  powerful  remote 
sensing  technique  for  characterizing  the  microphysics  of  contrail  particles.  Since  young 
contrails  often  consist  of  relatively  small  ice  crystals,  the  quantitative  interpretation  of  lidar 
measurements  requires  accurate  theoretical  computations  of  6  for  polydisperse,  randomly 
oriented  nonspherical  particles  with  size  parameters  ranging  from  zero  to  at  least  several 
tens,  thus  ruling  out  most  of  the  currently  available  numerical  techniques.  We  used  the 
recently  improved  T-matrix  method  and  computed  6  for  polydispersions  of  randomly 
oriented  ice  spheroids,  circular  cylinders,  and  Chebyshev  particles  with  sizes  typical  of 
young  contrails.  We  showed  that  ice  crystals  with  effective  radii  as  small  as  several  tenths 
of  a  micron  can  already  produce  6  exceeding  0.5  at  visible  wavelengths.  This  may  explain 
the  frequent  occurrence  of  large  6  values  for  very  young  contrails.  We  also  showed  that 
observed  increases  of  6  with  the  contrail's  age  can  be  explained  either  by  a  rapid  increase 
of  the  particle  size  parameter  from  essentially  zero  to  about  5  Or  by  assuming  that  the 
contrail  particles  originate  as  perfect  spheres  and  then  acquire  a  certain  degree  of 
asphericity. 

This  research  was  summarized  in  Ref.  4. 
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nonspherical  particles.  Second,  we  use  exact  T-matrix  computations  and  the  Kirchhoff  approximation  to  show  that  the  6- 
function  transmission  peak  predicted  by  the  GO  approximation  for  hexagonal  ice  crysals  is  an  artifact  of  GO  completely 
ignoring  physical  optics  effects  and  must  be  convolved  with  the  Fraunhofer  pattern,  thereby  producing  a  phase  function 
component  with  an  angular  profile  similar  to  the  standard  diffraction  component.  Third,  we  have  used  the  improved  T-matrix 
method  to  compute  the  linear  depolarization  ration  for  polydispersions  of  randomly  oriented  ice  spheroids,  circular  cylinders, 
and  Chebyshev  particles  writh  sizes  typical  of  young  contrails.  We  have  shown  that  ice  crystals  with  effective  radii  as  small 
as  several  tenths  of  a  micron  can  already  produce  5  exceeding  0.5  at  visible  wavelengths. 
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